Mon. Not. R. Astron. Soc. OOO.[THT2lC2002') Printed 27 May 2010 (MN L«?eX style file v2.2) 



Shear-driven and diffusive helicity fluxes in aQ dynamos 

G. Guerrero^* and P. Chatterjee^ and A. Brandenburg^'^ 

^ Nordita, AlbaNova University Center, Roslagstullsbacken 23, SE 10691 Stockholm Sweden 

^Department of Astronomy, AlbaNova University Center, Stockholm University, SE 10691 Stockholm, Sweden 



ABSTRACT 

We present nonlinear mean-field aO, dynamo simulations in spherical geometry with 
simplified profiles of kinematic a effect and shear. We take magnetic helicity evolution 
into account by solving a dynamical equation for the magnetic a effect. This gives 
a consistent description of the quenching mechanism in mean-field dynamo models. 
The main goal of this work is to explore the effects of this quenching mechanism 
in solar-like geometry, and in particular to investigate the role of magnetic helicity 
fluxes, specifically diffusive and Vishniac-Cho (VC) fiuxes, at large magnetic Reynolds 
numbers [R-m]- For models with negative radial shear or positive latitudinal shear, the 
magnetic a effect has predominantly negative (positive) sign in the northern (southern) 
hemisphere. In the absence of fluxes, we find that the magnetic energy follows an R^ 
dependence, as found in previous works. This catastrophic quenching is alleviated in 
models with diffusive magnetic helicity fluxes resulting in magnetic fields comparable 
to the equipartition value even for R^ — 10^. On the other hand, models with a 
shear-driven Vishniac-Cho flux show an increase of the amplitude of the magnetic 
field with respect to models without fiuxes, but only for R^ < 10**. This is mainly a 
consequence of assuming a vacuum outside the Sun which cannot support a significant 
VC flux across the boundary. However, in contrast with the diffusive fiux, the VC flux 
modifies the distribution of the magnetic field. In addition, if an ill-determined scaling 
factor in the expression for the VC fiux is large enough, subcritical dynamo action is 
possible that is driven by the action of shear and the divergence of current helicity 
flux. 

Key words: magnetic fields — MHD — hydrodynamics - turbulence 
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, 1 INTRODUCTION in the nonlinear phase of the dynamo process, conservation 

. . , T r of magnetic helicity gives rise to a magnetic a effect (om) 

>A crucial point m the study or astrophysical dynamos is .,, . -^ ^ ^i • i x- j_ ■■, ^- i ^ ^i 

-^ , , , ^T with a sign opposite to the inductive contribution due to the 

•i-H ,. to understand the mechanism by which they saturate. Nev- i i- i .• . .i i . ,■ rr , a .i i 

k> r 1 , helical motions, i.e., the kinematic a effect. As the produc- 

J\. ' ertheless, a consistent description or this process has rarely ,. „ , , 7-,,if.ii r.i 

^ , ^ . , . , , . n , r , /^ rr,T^\ , ,. , tion or om depends on Km, the fanal value ot the magnetic 

^ been considered in mean-held dynamo (MtD) modeling and j:: i i i i i i r n j_i i j tt i 

jrt ■ A held should also follow the same dependence. However, real 

. - - I only a heuristic description is often used. An important phe- ii-ii^i- j_ij ^ i^^i ^ 

■' , , , , , astrophysical bodies are not closed systems, but they have 

nomenon happens when the dynamo operates m closed or , i • ^i ^ n a c _l- i i- -^ 

J i- open boundaries that may allow a nux ot magnetic helicity. 

periodic domains: the turbulent contribution to the dynamo mi i i i- r i- i i- -^ -i- ^ _li 

-^ ihe shedding ot magnetic helicity may mitigate the catas- 

equation, i.e., the a effect, decreases for large values of the , i • i . 

f trophic Q quenching, 
magnetic Reynolds number. This process is known as catas- 
trophic quenching and can pose a problem in explaining the These ideas have been tested in direct numerical 

generation of magnetic field in late type stars hke the Sun or simulations (DNS) in both local C artesian and global 

the Galaxy, where _Rm could be of the order of 10^ or 10^^, spherical domains. In the former (IB randenburj l2005l : 

respectively. iKapvla. Korpi and B randen burg |2008| ) it has been clearly 

In the last few years the nature of the catas- shown that open boundaries (e.g. vertical field boundary 

trophic quenching has been identified as a consequence conditions) lead to a faster saturation of a large-scale mag- 

of magnetic hehcity conservation (for a review see netic field compared with cases in closed domains (per- 

iBrandenburg fc Subramanianl2005al ). It has been found that feet conductor or triple-periodic boundary conditions). In 

the latter, it has been found that it is possible to build 

up large-scale magnetic fields either w ith forced turbulence 

* E-mail: guerrero@nordita.org (GG) l|Brandenburgjl2005l : lMitra et al.ll2010H ) or with convectively 

© 2002 RAS 



2 G. Guerrero, P. Chatter jee and and A. Brandenburg 



drive n turbulence (e.g., I Brown et al.l l2010l : iKapvla et al.l 
I2OI0I ). These models generally used vertical field boundary 
conditions. 



In flux- transport dynamos (IDikpati &: Charbonneaul 

I1999I : iGuerrero fc de Gouveia Dal Find l2008l ) 

as well as in interface dynamos of t he so- 
lar cycle (e.g. MacGregor fc Charbonneaul 119971 : 
ICharbonneau fc MacGregoiJ 119971 ) the quenching mech- 
anism has been considered either through an ad hoc 
algebraic equation or b y phcnomcnological considerations 
l|Chatteriee. Nandv fc C houdhuri 20dl), but most of the 
time the models do not consider the effects of magnetic 
helicity conservation. An exception i s the r ecent paper by 
IChatteriee. Brandenburg fc Guerrero! (|2010l '). where these 
effects have been considered in the context of an interface 
dynamo. 

In general the magnetic helicity depends on time, so 
it is necessary to solve an additional dynamical equation 
for the contribution of the small-scale field to the magnetic 
helicity together with the induction equation for the mag- 
netic field. In the past few years, some effort has already 
been made to consider this dynamical saturation mecha- 
nism in M FD models like in the ID a^ dynamo models pre - 
sented in [Brandenburg. Candelaresi fc Chatteried l|2009t ). 
in axisymmetric model s in cylindrical geome try for the 
galactic afi dynamo (jShukurov et al.l I2OO6I ), and also 
in models with spheric al geometry for an o? dynamo 
IjBrandenburg et al.ll2007l ') . The role of various kinds of mag- 
netic helicity fluxes have be en explored in severa l papers 
llBraJidenburg. CandelaresTfc Chatterieeii2009i : .Zhang et al.l 
l2006l : IShukurov et al.ll2006l ). 



Our ultimate goal is to develop a self-consistent MFD 
model of the solar dynamo, with observed velocity profiles 
and turbulent dynamo coefficients computed from the DNS. 
This is a task that requires intensive efforts. Hence we shall 
proceed step by step, starting with simple models and then 
including more realistic physics on the way. In this work we 
will study the effects of magnetic helicity conservation in 
simplified aO. dynamo models for a considerable number of 
cases. More importantly, we shall perform our calculations in 
spherical geometry, which is appropriate for describing stel- 
lar dynamos, with suitable boundary conditions, and con- 
sidering shear profiles which are a simplified version of the 
observed solar differential rotation. We shall also explore 
how magnetic helicity fluxes affect the properties of the so- 
lution. Two classes of fluxes are considered in this paper: a 
diffusive flux and a shear-driven or Vishniac-Cho (hereafter 
VC) flux l|Vishniac fc Choi 120011 ). We consider models with 
either radial or latitudinal shear. The effects of meridional 
circ ulation will be investigated in detail in a com panion pa- 
per (|Chatteriee. Guerrero fc Brandenburell2010l '). 

This paper is organized as follows: in Section [5] we de- 
scribe the basic mathematical formalism of the aO, dynamo, 
give the formulation of the equation for qm and also justify 
the fluxes included. In Section [3] we describe the numerical 
method and then, we present our results in SectionUJstarting 
from a dynamo model with algebraic quenching to models 
with dynamical a quenching and different fluxes. Finally, we 
provide a summary of this work in Section [Sj 



2 THE aQ. DYNAMO MODEL 

In mean-field dynamo theory, the evolution of the magnetic 
field is described by the mean-field induction equation, 

^ = Vx (C7xB + f -r,„,VxB), (1) 

at 

where B and U represent the mean magnetic and ve- 
locity flelds, respectively, ?7m is the molecular diffusivity, 
£ = aB — tji^qJ is the mean electromotive force obtained 
using a closure theory like the flrst order smoothing approx- 
imation, where £ gives the contribution of the small-scale 
components on the large-scale field, a is the non-diffusive 
contribution of the turbulence, r]t is the turbulent magnetic 
diffusivity, J = V x B/no is the mean current density, and 
Ho is the vacuum permeability. 

In spherical coordinates and under the assumption of 
axisymmetry, it is possible to split the magnetic and the ve- 
locity fields into their azimuthal and poloidal components, 
B = Be^ -I- V X {Ae.^) and U = rsm6Q,e.^ + tip, respec- 
tively. For the sake of simplicity we shall not consider the 
meridional component of the flow, i.e. Up = 0. Then, the 
toroidal and poloidal components of equation (|l| rnay be 
written as 

^ = sBp ■ Vn - [V77 X (V X Be^)]^ + r^D^B , , (2) 

^ = aB + riD''A, 
at 

where D^ 



(3) 



V — s~ is the diffusion operator, rj — rim + rit, 
s — r sin 9 is the distance from the axis, and Bp = V x (Ae^) 
is the poloidal fleld. 

The two source terms in equations (|2]| and (|3]l, sBp- Vf2 
and aB, express the inductive effects of shear and turbu- 
lence, respectively. The relative importance of these two ef- 
fects may be quantified through the non-dimensional dy- 
namo numbers: Cn — AQ,L lr\\ and Ca = aoL/rit, where 
Ar2 is the angular velocity different between top and bot- 
tom of the domain. Note that equations @ and (|3I are valid 
only in the limit Cn ^ Ca, known as aQ. dynamo. 

The inductive effects of the shear may be understood as 
the stretching of the magnetic field lines due to the change 
in the angular velocity between two adjacent points. On the 
other hand, the kinematic a-effect is the consequence of heli- 
cal motions of the plasma which produce screw-like motions 
in the rising blobs of the magnetic field. Using the first order 
smoothing approximation it may be expressed as: 



Qk 



(4) 



where, r is the correlation time of the turbulent motions 
and oj = V X n is the small-scale vorticity. The saturation 
value of the magnetic field may be obtained by multiply- 
ing Ok by the quenching function /q = (l + B^/Bcq) , 
which saturates the exponential growth of the magnetic field 
at values close to the equipartition field strength given by 
_Beq = (/iopu^). This form of algebraic quenching was intro- 
duced heuristically (see, e.g. Stix 1972) and has been often 
used as the standard quenching mechanism in many dynamo 
simulations. However, it does not give information about the 
back reaction process and is independent of any parameter of 
the system like the magnetic Reynolds number. A consistent 
description of the quenching mechanism will be presented in 
the following section. 
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2.1 Dynamical a effect 

Recently, it has been demonstrated that when the amplitude 
of the magnetic field reaches values near the equipartition, 
the Q-efi'ect is modified by a magnetic contribution, the so 
called magnetic a efi'ect, denoted by ajvi- It is usually the 
case that au has a sign opposite to ok resulting thus in the 
satura tion of the magnetic field. IPouauet. Frisch fc Leorad 
l| 19761 ) have shown that am is proportional to the small- 
scale current helicity of the system, hence it is possible to 
write a as a sum of two contributions, one from the fluid 
turbulence and other from the magnetic field, as follows: 



"K + OM = —^TUJ ■ u + \tj ■ b/p , 



(5) 



where p is the mean density of the medium, assumed 
here as a constant, and J = V x b//io is the current 
density of the fiuctuating field. The mathematical expres- 
sion that describes the evolution of om rnay be obtained 
by taking into account the m agnetic helicity evolution 
l|Blackman fc Brandenburej|2002l ). which leads to: 



dam „ , 2 



g B om 
R2 R 



V -Tc 



(6) 



where fcf = 2it/{L — Vc) with re = O.TLq is a suitable 
choice for the wave number of the forcing scale, the mag- 
netic Reynolds number _Rm = Vt/Vtn and J^a is the fiux of 
the magnetic a effect related to the fiux of the small-scale 
magnetic helicity, Ff through: 



-'--'Gq 



(7) 



According to previous authors ajvi has a finite value in the 
interior of the domain in absence of fiuxes {J^a = 0), and 
its sign is usually opposite to the sign of qk in such a way 
that the final amplitude of the total a-effect decreases, and 
so does the final value of the magnetic energy. 



2.2 Magnetic helicity fluxes 

Recently it has been pointed out that the catastrophic 
quenching could be alleviated by allowing the fiux of small- 
scale magnetic (or current) helicity out of the domain, so 
that the total magnetic helicity inside need not be con- 
served any longer. Alternately, we may introduce those 
fiuxes in the equation for am', see equation ([7|. Several can- 
didat es have been proposed for the helicity fiuxes in the 
past jKleeorin fc Roga chcvskii' 1999; 'Vish niac fc Chdl200ll : 
ISubramanian fc Brandenburg 2004 ). Amongst them are the 
fiux of magnetic helicity across the iso-rotation contours, ad- 
vective and diffusive fiuxes and also the explicit removal of 
magnetic helicity in processes like coronal mass ejections or 
galactic fountain flows, for the case of the galactic dynamo. 

From the mathematical point of view, the nature of the 
fiux terms in the equation for qm has not been demonstrated 
with sufficient rigor. However, several DNS have pointed to 
its existence. 

Firstly, the shearing box c onvec tion simulations of 
iKapvla. Korpi and Brandenburd (|200a) showed that in the 
presence of open boundaries, the large-scale magnetic field 
grows on temporal scales much shorter than the dissipa- 
tive time scale. They concluded from this that open bound- 
aries may allow the magnetic helicity to escape out of the 



system. These experiments seem to be compatible with the 
fiux proposed by [Vishniac fc Cho (200l|), whose functional 
form may be expressed as fsee | Subramanian fc Brandenburel 
I2OO4I : iBrandenburg fc SubramanianI l2005bl . for further de- 
tails): 



Ti — Cyc>^ijiStkBjBk , 



(8) 

where Sik = ^{Ui^k + Uk,i) is the mean rate of strain tensor 
and Cvc is a non-dimensional scaling factor. As we assume 
Up = 0, this fiux has the following three components: 






^d, 



Cvc [S^rBgB^ + Se</)(Be — B^)\ 
Cvc [S^gBrBg + Sr4>{B^ — B^ 
Cvc pe^Br-B^ — Sr,j,{Bg — B^)\ 



(9) 
(10) 

(11) 



360 



with S</,r = Sr0 — rsai6{d^/dr)/2 and Sg 

sme{dQ/de)/2 

Secondly, iMitra et al.l (|2010al ) performed a^ dynamo 
simulations driven by forced turbulence in a box with an 
equator. They found that the diffusive fiux of om across the 
equator can be fitted to a Fickian diffusion law given by. 



Tt> = 



,(r)VQM 



(12) 



They also computed the numerical value of this diffusion co- 
efficient, and found it to be of the order of turbulent diffusion 
coefficient. They also found that the time averaged fiux is 
gauge independent. Both results were later corroborated by 
simulations without equator, but with a decline of kinetic 
helici ty toward the boundaries (jHubbard fc Brandenburd 
I2OIOI ). 

Additionally, magnetic helicity may be advected by the 
mean velocity with a fiux given by T^^ = omC/, or it may 
be expelled from the solar interior by coronal mass ejections 
(CMEs) or by the solar wind. This fiux, J-cme, rnay account 
for ~ 10% of the total helicity generated by the solar differ - 
ential rotation, as estimated bv lBerger fc RuzmaikinI (|2000|). 
It can be modeled by artificially removing a small amount 
ofoM every r time (IBrandenburg. Candelaresi fc Chatteried 
120091 ) ■ or also by a radial velocity field that mimics the solar 
wind. 

The total flux of magnetic helicity may be written as 
the sum of these contributions, 



T = T^ 



Ty)+T^ 



(13) 



Since in this dynamo model we do not include any compo- 
nent of the velocity field other than the differential rotation, 
in this study we will consider only the first two terms on the 
rhs of equation (|13|) . 



3 THE MODEL 

We solve equation ((2]), ((3]) and ([6]) for A, B and qm in the 
meridional plane in the range 0.61/ ^ r ^ L and ^ 6 ^ tt. 
We consider two different layers inside the spherical shell. In 
the inner one the dynamo production terms are zero and go 
smoothly to a finite value in the external layer. The mag- 
netic diffusivity changes from a molecular to a turbulent 
value from the bottom to the top of the domain. This is 
achieved by considering error function profiles for the mag- 
netic diffusivity, the differential rotation, and the kinetic a 
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effect, respectively (see Fig. [T]): 

aK{r,9) = Ca { —] 0{r,ri,wi)cos( 
\VtJ 



(14) 
(15) 

(16) 



where 0{r,ri^2,'w) = i [1 + erf {(r — r-i,2)/'"^i}], with ri = 
0.7Lo, ra = 0.72Lo and wi = 0.025Lo. We fix Cn = -10^ 
and vary Ca- 

The boundary conditions are chosen as follows: at the 
poles, 9 = 0,TT, we impose A = B — 0; at the base of the 
domain, we impose a perfect conductor boundary condition, 
i.e. A = d{rB)/dr = 0. Unless noted otherwise, we use at 
the top a vacuum condition by coupling the magnetic field 
inside with an external potential field, i.e., (V^ — s~^)A = 0. 
A good description of the numerical i mplementation of this 
boundary condition may be found in lDikpati fc Choudhuril 
(| 19941 ). 

The equations for A and B axe solved using a second- 
order Lax-Wendroff scheme for the first derivatives, and cen- 
tered finite differences for the second-order derivatives. The 
temporal evolution is computed by using a modifie d ver - 
sion of the ADI method of Pcaccman & Rachfordj ()l955l l 
as explained in iDikpati fc Charbonneaul (| 19991 ). This nu- 
merical scheme has been used previously in several works 
on the flux-transport dynamo and the results were found 
to be in good agreement with those using other numerical 
techniques ([Guerrero fc de Gouvcia Dal Pino 2007], 120081 : 
iGuerrero. Dikpati fc de Gouveia Dal Pindl2009l ). 

In the absence of magnetic helicity fluxes, equation ([6]) 
for QfM corresponds to an initial value problem that can be 
computed explicitly. However, as we are going to include a 
diffusive flux, we use for am the same numerical technique 
used for A and B. All the source terms on the right hand 
side of equation ^ are computed explicitly. We have tested 
the convergence of the solution for 64^, 128^, and 256^ grid 
points. For cases with small Rm, there are no signiflcant 
differences between different resolutions, but for high i?m, 
&4? grid points is insufficient to properly resolve the sharp 
diffusivity gradient. A resolution of 128^ grid points is a 
good compromise between accuracy and speed. 



4 RESULTS 

4.1 aO. dynamos with algebraic quenching 

In order to characterize our afi dynamo model we start by 
exploring the properties of the system when the saturation is 
controlled by algebraic quenching with fq = (l + _B^/_B^q)~^- 
We found that, with the profiles given by equations (I14II - 
()16p. Fig. m the critical dynamo number is around 2 x 10* 
(i.e., Cq = 1.975). The solution for the model is a dynamo 
wave traveling towards the equator since it obeys the Parker- 
Yoshimura sign rule (see Fig.[2|. In this case, the maximum 
amplitude of the magnetic field depends only on the dy- 
namo number of the system, CaCn, as can be seen in the 
bifurcation diagram in Fig.O The quenching formula is here 
independent of Rm, so the saturation amplitude is also in- 
dependent on Rm- 
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Figure 1. Profiles of the dynamo ingredients, cjk (solid line), 
dQ/dr (dashed line) and rjt (dot-dashed line). All the values are 
normalized to their maximum value. 
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Figure 2. Time-latitude butterfly diagram for the toroidal com- 
ponent of the magnetic field at r = 0.72Lo, for an aO, dynamo 
model with Cn = —10'' and Ca = 2.5. 



4.2 aO. dynamos with dynamical quenching 

In this section we consider dynamo saturation through the 
dynamical equation for qm described in Section [2]T] In this 
models we distinguish three different stages in the time evo- 
lution of the magnetic field: a growing phase, a saturation 
phase and a final relaxation stage (see panels a, b and c of 
Fig. O . The magnetic field is amplified from its initial value, 
5 X lO^^Boq, following an exponential growth. From the ear- 
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Table 1. Summary of main parameters and results of the numerical simulations. 



Run 


Cc 


^m 


Ka (r)t) 


Cvc 


-t>rms/-Dcq 


T{Ll/ri,) 


t (Ll/vt) 


CaC 


1.975 


10 


- 


- 


0.0008 


0.0486 


1.0 


Ca2.0 


2.0 


10 


- 


- 


0.15 


0.0484 


1.0 


Ca2.1 


2.1 


10 


- 


- 


0.33 


0.0477 


3.0 


Ca2.2 


2.2 


10 


- 


- 


0.45 


0.0471 


3.0 


Ca2.3 


2.3 


10 


- 


- 


0.56 


0.0464 


3.0 


Ca2.4 


2.4 


10 


- 


- 


0.64 


0.0460 


3.0 


Ca2.5 


2.5 


10 


- 


- 


0.69 


0.0455 


3.0 


RmlO 


2.5 


10 


- 


- 


0.21 


0.0422 


4.0 


Rm50 


2.5 


50 


- 


- 


0.25 


0.0446 


4.0 


Rmle2 


2.5 


100 


- 


- 


0.2 


0.0455 


4.0 


RmleS 


2.5 


103 


- 


- 


0.07 


0.0464 


4.0 


Rm2e3 


2.5 


2x10^ 


- 


- 


0.05 


0.0464 


6.0 


Rm5e3 


2.5 


5x10^ 


- 


- 


0.03 


0.0468 


15.0 


Rmle4 


2.5 


10* 


- 


- 


0.02 


0.048 


15.0 


DRmlO 


2.5 


10 


0.005 


- 


0.22 


0.0422 


4.0 


DRmSO 


2.5 


50 


0.005 


- 


0.26 


0.0446 


4.0 


DRmle2 


2.5 


100 


0.005 


- 


0.20 


0.0455 


4.0 


DRmle3 


2.5 


10^ 


0.005 


- 


0.09 


0.0460 


4.0 


DRmle4 


2.5 


10* 


0.005 


- 


0.06 


0.0457 


5.0 


DRmle5 


2.5 


10^ 


0.005 


- 


0.05 


0.0460 


7.0 


DRmle6 


2.5 


lO^ 


0.005 


- 


0.05 


0.0457 


8.0 


DRmle7a 


2.5 


lO'' 


0.001 


- 


0.026 


0.0457 


20.0 


DRmleTb 


2.5 


10^ 


0.005 


- 


0.05 


0.0460 


10.0 


DRmle7c 


2.5 


10^ 


0.01 


- 


0.073 


0.0460 


10.0 


DRmleTd 


2.5 


10^ 


0.03 


- 


0.12 


0.0460 


8.0 


DRmleTe 


2.5 


10^ 


0.05 


- 


0.15 


0.0457 


4.0 


DRmle7f 


2.5 


lO'' 


0.1 


- 


0.20 


0.0460 


4.0 


DRmle7g 


2.5 


10^ 


1.0 


- 


0.54 


0.0458 


4.0 


DRmlerh 


2.5 


10^ 


5.0 


- 


1.23 


0.060 


4.0 


DRmle7i 


2.5 


10^ 


10.0 


- 


1.76 


0.0457 


4.0 


VCa 


2.5 


10^ 


- 


0.002 


0.032 


0.0449 


4.0 


VCb 


2.5 


103 


- 


0.01 


0.02 


0.0442 


4.0 


VCc 


2.5 


10^ 


- 


-0.002 


0.02 


0.0447 


4.0 


VCd 


2.5 


104 


- 


-0.002 


- 


- 


4.0 


VCD 


2.5 


103 


0.1 


0.001 


0.11 


0.0446 


4.0 


ReleSe 


2.5 


10^ 
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0.023 


0.0282 
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0.004 


0.04 


0.033 
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0.004 
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Figure 3. Magnetic field average amplitude as a function of Ca 
using an algebraic quenching function that is independent of iJm ■ 



liest stages of the evolution we notice the growth of cxyi 
with values that are predominantly negative in the northern 
hemisphere and positive in the southern hemisphere. The 
latitudinal distribution of qm is fairly uniform in the ac- 
tive dynamo region, spanning from the equator to ~ 60° 
latitude. The radial distribution exhibits two narrow layers 
where the sign of om is opposite to the dominant one devel- 
oping at each hemisphere. These are located at the base of 
the dynamo region (r ~ O.TLo) and at a thin layer near to 
the surface (r > 0.95). In the equation for the magnetic a 
effect, equation (|6}, the production term is proportional to 
£ ■ B = aB — r^tfioJ ■ B. The first component of this term 
has the same sign as Qk, which in general is positive in the 
northern and negative in the southern part of the domain. 
The minus sign in front of the right hand side of equation 
(O defines then the sign of om- However, at the base and at 
the top of the dynamo region, OfK — >■ and B — >■ 0, respec- 
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tively. The term rjtJB is the only source of a-M and leads 
to the formation of these two thin layers. 

The space-time evolution of au depends on the value 
of the magnetic Reynolds number. For small i?m, the decay 
term in equation [6] (i.e. the second term in the parenthesis) 
becomes important, so that there is a competition between 
the production and decay terms resulting in an oscillatory 
behavior in the amplitude of the magnetic a effect, as is 
indicated by the vertical bars in the middle panel of Fig. [T] 
The period of these oscillations is the half the period of the 
magnetic cycle. With increasing Rm, the amplitude of the 
oscillations decreases such that for Rm ^ 10'^, om is almost 
steady. 

The morphology of the magnetic field corresponds to 
a multi-lobed pattern of alternating polarity (left panels of 
Fig. [5]). These lobes are radially distributed in the whole 
dynamo region with maximum amplitude at the base of this 
layer. The poloidal magnetic field follows a similar pattern 
with lines that are open at the top of the domain due to the 
potential field boundary condition. There is a phase shift be- 
tween toroidal and poloidal components which we have esti- 
mated to be ^ 0.4-K. The model preserves the initial dipolar 
parity during the entire evolution. 

The evolution of qm traces the growth of the magnetic 
field, but its final value depends on the magnetic Reynolds 
number. For small Rm, after saturation, om reaches a steady 
state, but for large Rm, its relaxation is modulated by over- 
damped oscillations. The relaxation time is proportional to 
Rm, which means that for Rm S> 1 the simulation must run 
for many diffusion times. The differences in the relaxation 
time observed for qm reflects the evolution of the magnetic 
field, as is shown in Fig. ID 

We observe that the rms value of the magnetic field 
remains steady during the saturation phase for Rm < 10^. 
For 10^ < Rm < 10'^ , a bump appears in the curve of mag- 
netic field evolution, followed by the relaxation to a steady 
value, whereas for Rm > 10^, the magnetic energy shows 
over-damped relaxations with a final energy proportional to 
Rm~ ^ as has been previously reported ( Brandenbu rg et al.l 
120071 ). These oscillations in the time evolution plot of 
the averaged magnetic field have been reported in mean 
field dynamo simulations including the dynamical a-effect 
([Brandenburg & Subramanian 2005b). 

Not many DNS of aQ, dynamo exist so far in the lit- 
erature with Rm ^ 100 in order to compare with our 
results. However, in the local aQ dyn amo simulations of 
iKapvla, Korpi and Brandenburg l|2008h . a rapid decay of 
the magnetic field seems to occur after the initial satura- 
tion for moderate values of Rm- This decay forms a bump in 
the curve of the averaged magnetic field (see their Fig. 14), 
similar to the bump that we obtain for lO'^ < Rm < 10^. 

For reasons of clarity in the Fig. 0] we do not show the 
entire time evolution of each simulation with Rm > 10'^. 
The total evolution time as well as the final value of the 
magnetic field of each simulation are shown in the Table [T] 
For magnetic Reynolds numbers above 2 x 10'*, the initial 
kinematic phase is followed by a decay phase during which 
the total a effect goes through subcritical values and then 
the dynamo fails to start again. 

In Fig. [5] we present the meridional distribution of the 
magnetic field (left panel), um (middle panel) and the to- 
tal a (right panel), in normalized units, for the three dif- 
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Figure 4. Time evolution of the averaged mean magnetic field for 
different values of Rm- Note that for Rm > 10'^, we have allowed 
the simulations to evolve more than 4 diffusion times, as indicated 
in Table [1] 



ferent stages of evolution corresponding to the early kine- 
matic phase, the late kinematic phase and the saturated 
phase. These snapshots correspond to the simulation with 
Rm = 10^ (Run Rmle3 in Table [1]). The multi-lobed pat- 
tern of the toroidal field represented with filled contours re- 
mains unchanged during the evolution even though its am- 
plitude increases. The same occurs for the poloidal compo- 
nent, shown by continuous and dashed streamlines for posi- 
tive and negative values, respectively. 

The magnetic a effect (middle panels) is formed first at 
latitudes between ±30° and then it amplifies and expands 
to latitudes up to ~ ±60°. This makes the total a effect, 
initially similar to ok (Fig. [l] and top panel of Fig. [5^), 
smaller at lower latitudes in the central area of the dynamo 
region. At the bottom and at the top of the domain am and 
qk have the same sign making the total a larger. However, 
the global effect is a decrease of the dynamo efficiency. 



4.3 Diffusive flux for om 

In this section we consider a Fickian diffusion term in equa- 
tion (|12|) for qm. We consider a diffusion coefficient varying 
from 5 X 10~^rit to lOr^t in the dynamo region and with 
Ka = r/m in the bottom layer. In these cases, the initial evo- 
lution of um is similar to the cases presented in the previous 
section: negative (positive) values for um in the northern 
(southern) hemisphere, with narrow regions of opposite val- 
ues nearby the regions where qk = or B = 0. However, at 
the later stages, um is much more diffuse in the entire do- 
main and has only one sign in each hemisphere. This is the 
result of cancellation of om with opposite signs occurring 
in each hemisphere due to radial diffusion. Contrary to the 
cases without fluxes, we now obtain finite values of _Bsat for 
large values of Rm, as can be seen in Fig.[Sl AH the cases de- 
picted in this figure correspond to Ka — O.OOSr^t. We notice 
that the final value of the magnetic field still remains small 
compared to the equipartition (^ O.lBcq), but it is clear 
that even this very modest diffusion prevents the a effect 
from being catastrophically quenched. This is also evident 
from the top panel of Fig.[7l where we plot the final strength 
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a) Magnetic field 





Total a 




-0.01-0.01 0.00 0.01 0.01 -0.00-0.00 0.00 0.00 0.00 -2.60-1.25 0.00 1.25 2.50 





-0,21-0.11 0.00 0.11 0,21 -2,50-1,25 0,00 1.25 2.50 




-0.11-0.06 0.00 0.06 0.11 -0.84-0.42 0.00 0.42 0.84 -2.50-1.25 0.00 1.25 2.50 



Figure 5. Meridional snapshots of three different phases of evolu- 
tion of the dynamo model with dynamical quenching, a) i = 0.25 
(L2/r,t), b) t = 0.5 (Llf-nt) and c) t = 2.0 {Lg/r?t). The left 
panel shows the contours of toroidal magnetic field in color scale, 
and positive (negative) poloidal magnetic field lines in continu- 
ous (dashed) lines. The central panel shows the distribution of 
oM) and the right panel shows the distribution of the residual a. 
All values are in non-dimensional units (i.e., B/Bcq), so that the 
color scale is different for each figure as indicated in the respective 
color bar. 



dissipative flux. In the middle and bottom panels of the Fig. 
[7] we compare the behavior of the normalized om, at a given 
point inside the dynamo region, and also the time period, T, 
of the dynamo for models with and without fluxes. In both 
panels it is clear that for Rm above ~ 10^, am and T reach 
a saturated value. 

Besides its dependence on Rm, the evolution of au de- 
pends also on Ka ■ For models with k^ <ti rjt, the evolution of 
qm relies on Rm, but for Ka ^ 0.l77t, the dissipation time of 
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Figure 6. The same that Fig. |4] but for simulations including 
a diffusive flux of ayi- AH the simulations correspond to Ka = 
0.005»7t. 



um becomes comparable to, or even shorter, than the period 
of the dynamo cycle. This results in a-M becoming oscilla- 
tory, as shown in the bottom panel of Fig. [S] The amplitude 
and the period of these oscillations depend on the value of 

In the top panel of Fig. [8] we show the final value of 
the averaged mean magnetic field as a function of Kq. We 
observe that for Ha in the range (0.1-1) rjt, the value of Brms 
remains between 20% and 60% of the equipartition, a value 
similar to the one obtained in the simulations using alge- 
braic a quenching (Section 14. 11 Fig. (Sj. For Ka > rjt, super- 
equipartition values of the magnetic field may be reached. 
This is because larger values of Kc result in oscillations of 
Qm with larger amplitude, such qm niay locally change its 
sign, increasing the value of the total a in each hemisphere 
and thereby enhancing the dynamo action. Such high val- 
ues of the diffusion of the magnetic helicity are unlikely in 
nature. 



4.4 The Vishniac-Cho flux 

Our next step is to explo r e the magnetic helicity flux pro- 
posed by Vishniac fc Choi (|200ll ) in the form given by equa- 
tion ((8|. For the moment we set Ka = 0. In a previous study 
on the effe cts of the VC flux in a MFD model in C artesian co- 
ordinates, |Brandenbujg_&_SuhramaB^^ (|2005a ) found that 
there exist a critical value for the parameter Cvc above 
which there is a runaway growth of the magnetic fleld that 
can only be stopped using an additional algebraic quenching 
similar to the one used in Section |4Tl] They found that this 
critical value, Cvc*, diminishes with increasing the amount 
of shear. Since we have used a strong shear (Cn ~ —10^) we 
use nominal values of Cvc = 10~^, but without any alge- 
braic quenching. 

The term V • jFvc develops a multi-lobed pattern which 
travels in the same direction as the dynamo wave, this con- 
firms that the VC flux follows the lines of iso-rotation. From 
equation ((8]), we see that the VC flux is proportional to the 
magnetic energy density. In the present case, with Cn ^ Ca, 
the spatial distribution of V -J-vc/B^^ is dominated by the 
terms involving B^ in equations. I ^lll (this may be inferred 
from the left hand panels of Fig. llUb .). This results in a 
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Figure 7. Rm dependence of the averaged mean magnetic field 
(top), the temporal mean value of qm at r = O.8L01 S = 45° 
(middle) and the dynamo cycle period, T in diffusion time units 
(bottom). The continuous line present the result for simulations 
without «M diffusive flux (kq = 0) and the dashed line shows the 
results for k^ = 0.005»7t. The error lines in the middle panel in- 
dicate the maximum and minimum amplitudes in the oscillations 
of Qm S't that point. 



new distribution of qm, with concentrated regions of positive 
(negative) sign at low latitudes in the northern (southern) 
hemisphere, and a broad region of negative (positive) sign 
in latitudes between 20° and 60° latitude (see middle pan- 
els of Fig. I10|) . Surprisingly we find that the general effect 
of this flux is to decrease the final amplitude of the mag- 
netic field with respect to the case without any fluxes as can 
be seen in Fig. 1111 Note that we have until now used only 
the potential field boundary condition for the poloidal field. 



Figure 8. Top, final amplitude of the rms mean magnetic field 
for different values of k, in this case Rm = 10^- Bottom, final 
amplitude of om at r = O.8L0 and 6 = 45°. The error lines 
indicate the maximum and minimum amplitude in the oscillations 
of om at this point. 



Magnetic field 





Total a 




-0.75-0.37 0.00 0.37 0.75 -0,53-0.26 0,00 0,26 0,53 -2,50-1,25 0,00 1,25 2,50 

Figure 9. The same than Fig. \E\ but for a diffusive flux, with 
K = r]. The snapshot corresponds to t = 3.0 (Lo/'yt). 



When we consider both diffusive as well as VC fluxes, with 
Ka = 0.1?7t and Cyc = 10~^, we obtain a magnetic field of 
slightly larger amplitude compared to the case with only the 
diffusive flux (compare the value of firms in Runs DRmleS 
and VCD in Table [T} . However we may say from the but- 
terfly diagram of Fig. [T^] that the toroidal magnetic field 
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appears to be more concentrated at lower latitudes, where 
the sign of Qm is same as that of ax- 

With negative values of Cvc, it was found that the 
resulting profile of om is only weakly modified from cases 
without fluxes, though its value is reduced marginally such 
that the final amplitude of Brms is slightly larger. But even 
this contribution does not help in alleviating catastrophic 
quenching in models with large 7?m (see Fig. Ilip. 

Since VC fluxes transport helicity along lines of con- 
stant shear, it may be expected that they are more impor- 
tant in models with latitudinal shear, since in this case the 
magnetic helicity flux can travel either towards the bottom 
or the top boundaries, from where magnetic helicity can be 
expelled. For testing this possibility, we turn off the radial 
shear profile and consider a purely latitudinal solar-like dif- 
ferential rotation: 



n{r,e) = Cn 



Vt 



o r2 



eir,r2,wi){ns{e)~nc 



(17) 



where f2eq/27r — 460.7 nHz is the angular velocity at the 
equator, and fla (0) = ilcq + 02 cos^ + a^ cos'' gives the 
latitudinal profile, with a2/2n — —62.9 nHz and a4,/2TT = 
-67.13 nHz. 

In order for the dynamo to be slightly supercritical, 
as in the previous cases, we consider Cn = 5 x 10*. This 
dynamo solution corresponds to a dynamo wave produced 
at mid latitudes (~ 45°) that travels upwards (since Cn 
now is positive). As in the previous cases with radial shear, 
the distribution of V ■ J^vc/B^q is similar to that of the 
divergence of magnetic energy density (left hand panels of 
Fig. [To] b,c and d). If no fluxes are considered, the final 
amplitude of the mean magnetic field is ~ 0.03% of the 
equipartition value. In presence of VC fluxes, starting with 
Cvc = 10~^ for a model with Rm = 10^, we notice that 
the final magnetic field is twice as large as in the case with 
Cvc = 0. 

Our model becomes numerically unstable beyond 
Cvc = 10~^ due to appearance of concentrated regions of 
strong am- When VC and diffusive fluxes are considered si- 
multaneously, with Cvc ~ W~^ and Kc = 0.1?7t, the relaxed 
value of -Brms is only slightly below the value reached at 
the end of the kinematic phase (Fig. Illb ). In this case qm 
spreads out in the convection zone, as shown in Fig. llOfc . 
indicating that the effects of the VC flux are not important 
when compared with the diffusive flux. 

We repeated the calculation by considering the vertical 
field (VF) boundary condition, d{rBg)/d6 — 0, for the top 
boundary, instead of the potential field (PF) condition used 
throughout the rest of this work. Furthermore, in the mod- 
els with VF conditions the presence of the VC flux leads to 
an increase of -Baat by a factor of ~ 2 compared to the case 
without VC flux (see Fig. I lib ). It may be noted that qm 
shows regions of both positive and negative signs in each 
hemisphere (see Fig. IIOH ). Thus, the total a effect is in- 
creased locally to values well above the kinematic one. This 
implies that in the region around ±45° the dynamo action 
is driven by the magnetic a effect. A similar secondary dy- 
namo is found to be working for a different distribution o f 
shear and uk (|Chatteriee. Guerrero fc Brandenburell201(t ). 
As with PF boundary condition, large values of Cvc result 
in a numerical instability of the magnetic fleld in the simu- 
lation with VF. 
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Figure 11. Time evolution of the averaged mean magnetic field 
for different values of Cvc- ^) Radial shear, b) latitudinal shear 
with potential field boundary conditions and c) latitudinal shear 
with vertical field boundary conditions. The width of the differ- 
ent bands reflects the range over which the magnetic field varies 
during one cycle. Note that the cycle period is short compared 
with the resistive time scale on which the magnetic field reaches 
its final saturation. If not indicated, in all models i?m = 10^. The 
two dashed lines in the panel a) corresponds to Cvc = —0.002 for 
iim = 10^ and Rm = 10*. 



The main result of this section is that the VC flux does 
not alleviate catastrophic quenching of the dynamo for large 
values of _Rni (see the dashed lines in Fig. [TT]a and c). The 
reason for this may be related to the fact that the radial flux 
has components that are either proportional to Bg or to B^ 
(equation [5]). As B^ vanishes on the top boundary, and Bg 
is small, the VC flux is not able to dispose of qm across 
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a) Magnetic field ^-Fyc/B^ 



Total a 




-0,0-1-0,02 0,00 0,02 0,01 -0,50-0,25 0,00 0,25 0,50 -1,14-0,57 0,00 0,57 1,14 -2,52-1,26 0,00 1,26 2,52 

c) 




-0.11-0,05 0,00 0,05 0,11 -4,32-2,16 0,00 2,16 4,32 -3,72-1,86 0,00 1,86 3,72 -5,60-2,80 0,00 2,80 5,60 



Figure 10. Meridional snapshots of different models the in Table [T] a) model VCa, b) model VCga, c) model VCDg and d) model 
VCgvf. The contours (colors and lines) for the magnetic field have the same meaning than in Fig. [S] In this plot we have include a new 
column with the value of the VC component in the om equation, i.e., V ■ IFvc/I^cq- ^11 ^^^ snapshots corresponds to the relaxed state 
of evolution. (g 2002 RAS, MNRAS 000.ITHT21 
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a) Model Rmlc 
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Figure 12. Butterfly diagrams of toroidal field for runs without 
magnetic helicity flux (a) and with VC flux (b) for Rm = 10"' ■ 
Note the stronger concentration of magnetic field at lower lati- 
tudes in the presence of VC flux. 



the boundary. This might change if diffusive fluxes became 
important near the top or if a different boundary condition 
on B were apphed. 



5 CONCLUSIONS 

We have developed aQ dynamo models in spherical geome- 
try with relatively simple profiles of qk and shear {dil/dr 
and do,/ 86). We choose potential field (also vertical field 
in some cases) and perfect conductor boundary conditions 
for the top and bottom boundaries, respectively. We esti- 
mate the critical dynamo number by fixing Cn ~ —10'* and 
varying Ca while using algebraic quenching. 

Using a dynamo number, CnCa, that is slightly super- 
critical, we solve the induction equations for B and A to- 
gether with an equation for the dynamical evolution of the 
magnetic a effect or om- We find that for positive (negative) 
values of Ca in the northern (southern) hemisphere, om is 
mainly negative (positive), with narrow fractions of opposite 
sign in regions where ok or B are equal to zero. 

We find that the kinematic phase is independent of 
Rm- However for Rm > 10^ there exists a phase of 
relaxation post saturation in which the averaged mag- 
netic field oscillates about a certain mean. The larger 
the Rm, the more pronounced are the damped oscilla- 
tions and the longer is the relaxation time (Fig. |3|. The 



final value of the magnetic energy obeys a R^^ depen- 
dency (i?""'^ for magnetic field, Fig. [71), which is in 
agreement with earlier work ([Brandenburg fc SubramanianI 
[2005bl : [Brandenburg. Candelaresi fc Chatteriee|[2009l l. 

We argue that including equation ([SJ in MFD models 
is appropriate for describing the quenching of the magnetic 
field in the dynamo process. Since we observe large-scale 
magnetic fields at high magnetic Reynolds numbers in as- 
trophysical objects, there must exist a mechanism to prevent 
the magnetic field from catastrophic quenching. 

We have studied the role that diffusive and VC fluxes 
may play in this sense. Their contribution may be summa- 
rized as follows: 

(i) In the presence of diffusive fluxes, am has only one sign 
in each hemisphere (negative in the northern hemisphere 
and positive in southern) and is evenly distributed across 
the dynamo region (Fig. [9} . 

(ii) For Rm < 10^ the mean values of qm are similar to 
models without diffusive fluxes, whereas for Rm ^ 10^, ayi 
has smaller values that seem to be independent of Rm (see 
Fig. [3 middle). 

(iii) Even a very low diffusion coefficient, e.g. Kc, — 
O.OOlryt, causes Brms to depart from the Rm^'^ tendency and 
converge to a constant value which is then around 5% of 
the equipartition value for large values of Rm, but below the 
value of 10^ used in this study (dashed line in Fig. [7] top). 

(iv) Larger values of Ka result in larger final field 
strengths. 

(v) In models with only radial shear the Vishniac-Cho 
fiux contributes to qm with a component that travels in 
the same direction as the dynamo wave. This produces a 
different radial and latitudinal distribution of the magnetic 
a effect that also affects the distribution of the magnetic 
fields. However, it does not help in alleviating the quenching 
at high Rm- On the contrary, the larger the coefficient Cvc, 
the smaller is the resultant magnetic field. 

(vi) In models with only latitudinal shear the VC fiux 
travels radially outward but it remains concentrated at the 
center of the dynamo region. In a given hemisphere the re- 
sultant distribution of au has both positive and negative 
signs. The part of om that has the same sign as qk en- 
hances dynamo action. This effect is more evident in models 
with vertical field boundary conditions (Figs. fTOb -d'l. 

(vii) In models with vacuum and vertical field boundary 
conditions and Rm ~ 10^, the VC fiux increases the final 
value of the magnetic field by a factor of two compared to 
the case without any fluxes. 

(viii) The magnetic field in models with Rm ^ 10* and 
with non-zero VC fiux decays after the kinematic phase since 
the total a effect becomes subcritical (see dashed lines in 
Fig. [11] a and c) . 

(ix) Larger values of Cvc produce narrow bands of au 
which drives intense dynamo action in these regions. This 
positive feedback between the magnetic field and qm causes 
the simulation to become numerically unstable in the ab- 
sence of any other quenching effect. 

From the above results it is clear that diffusive fluxes are 
much more important in alleviating catastrophic quenching 
when compared to the Vishniac & Cho fluxes (in the form 
of equation [HJ for a large range of Rm- This is somehow 
intriguing since it is known from DNS that shear in do- 
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mains with open boundaries does indeed help in aUeviating 
the catastrophic quenching. It may be understood as a re- 
suh of the large value of Cn compared with Ca and also to 
the top boundary condition for the azimuthal magnetic fiel d 
l|Brandenburg||2005l : iKapvla. Korpi and BrandenburdbOOSD . 

The results presented above indicate that considerable 
work is still necessary in order to understand the role of 
larger-scale shear in transporting and shedding small-scale 
magnetic helicity from the domain. 

In snapshots of the meridional plane as well as in but- 
terfly diagrams we notice that the diffusive fluxes do not 
signiflcantly modify the morphology and the distribution 
of the magnetic fleld when compared with cases without 
fluxes or even with simulations with algebraic a quenching. 
On the other hand, for models with VC flux the distribu- 
tion of qm becomes different and so does the magnetic field. 
This is clear from the butterfly diagram shown in Fig. 112b . 
which exhibits a magnetic field confined to equatorial lati- 
tudes reminiscent of the observed butterfly diagram of the 
solar cycle. Even though this result corresponds to a simpli- 
fied model, it illustrates the importance of considering the 
dynamical a quenching mechanism for modeling the solar 
dynamo. Similar changes in the distribution of au and B 
are expected to happen when advection terms are included 
in the governing equations. 

In the simulations presented here, il and a effects 
are present in the same layers. An interesting ques- 
tion is whether the quenching of the dynamo is catas- 
trophic when both layers are segregated, as in the Parker's 
interface dynamo or the flux-transport dynamo mod- 
els. We addre ss this question in detail in two com- 
panion papers |Chatteriee. Brandenburg fc Guerrero! l20ld : 



IChatteriee. Guerrero fc BrandenburdboiOr ). 



We should notice that the back reaction of the 
magnetic fleld affects not only the a effect, but also 
the other dynamo coefficients, including the turbulent 
diffusivity. Contrary to quenching of a, the quenching 
of rjt may be considered through an algebraic quench- 
ing function ( see e.g. lYousef. Brandenburg fc Riidigen 



20031: 



Kapvla fc Brandenburg 



2009h . 



Guerrero, Dikpati fc de Gouveia Dal Pind (|2009l 'l have 



shown that in a flux-transport model these effects could 
affect properties of the models such as the final magnetic 
field strength and its distribution in radius and latitude. 
We leave the study of models with simultaneous dynamical 
a and 77 quenchings for a future paper. Solar-like profiles 
of differential rotation and meridional circulation along 
with dynamical a quenching will also be considered in a 
forthcoming paper. 
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